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A number of quantum algorithms have been performed on small quantum computers; these include Shor's 
prime factorization algorithm, error correction, Grover's search algorithm and a number of analog and 
digital quantum simulations. Because of the number of gates and qubits necessary, however, digital 
quantum particle simulations remain untested. A contributing factor to the system size required is the 
number of ancillary qubits needed to implement matrix exponentials of the potential operator. Here, we 
show that a set of tunneling problems may be investigated with no ancillary qubits and a cost of one 
single-qubit operator per time step for the potential evolution, eliminating at least half of the quantum gates 
required for the algorithm and more than that in the general case. Such simulations are within reach of 
current quantum computer architectures. 



Q 



uantum simulations on quantum computers are one of a set of algorithms that give exponential 
improvement in computational resources relative to the best classical algorithm 1 " 4 . 



Algorithms for studying many types of quantum field theory have been considered in both analog form in 
which a quantum Hamiltonian (typically manybody or multiple spin) is mapped either directly or via a suitable 
pulse-sequence to a computational Hamiltonian 5 " 21 and digital form in which a quantum system's Hamiltonian is 
split into free and interacting operators, then, using Trotter's formula, is simulated on a quantum computer 22 ' 23 . 

Small quantum simulations have already been realized on NMR 24 " 31 , atomic 32 ' 33 , ion trap 34 " 38 and photonic 39 ' 40 
quantum computers in both analog and digital forms. 

A concerted effort has also been made to investigate digital quantum particle simulations 3 for the simulation of 
chemical dynamics 41,42 . Algorithms for state preparation 43 , the simulation of temporal dynamics 41 and the mea- 
surement of observables 44 have all been developed. However, this type of simulation has remained untested due to 
the large number of gates and/or ancillary qubits needed to compute the kinetic and potential operators 41 ' 45 . In this 
paper, we focus on reducing the number of gates required for the simulation of temporal dynamics to the bare 
minimum, while still simulating interesting physics. 

For a review of the current state and outlook for quantum simulation, see 46 " 50 and 51 . 

The standard digital quantum simulation algorithm for a particle on a one-dimensional grid 3,45 encodes the 
position efficiently inn = log 2 N qubits, where AT is the size of the lattice of discretized particle locations, = kAx, 
k = 0, ... ,N — 1. The method uses a split operator approach to integrate a Schrodinger equation with a time- 
independent Hamiltonian that is first-order accurate in the time step, At 3 : 



\m)=e- m \^ init ) 



= e -i(V + K)t 



=( e -^ e -^< e ^))>,, it ). 

Higher order methods that give more accurate time integration have been developed 52 " 54 , but methods of order 
higher than two require more gates per time step. For simplicity, we will only consider first-order methods here. 
States in the qubit Hilbert space 

| l A>H---^^o>=---l^>IWh/'o>, 

where \i// t ) e {|0), |1)}, represent particle location in the binary representation, \x) = £/<=o,...,n - i The 
matrix exponential for the kinetic operator is calculated using a quantum Fourier transform (QFT) 
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Figure 1 | Particle probability distributions as a function of time for the 
first four steps of a two qubit simulation for v = 0 (free particle) and v = 
10 (particle in double well). The double- well potential and gray scale used 
to plot the probabilities are depicted to the right. 



--Fe 



where is a discrete Fourier transform operator and T is diagonal 
with entries proportional to —q 2 /2m, and q denotes the Fourier mode 
wavenumber. The resulting periodic, shift invariant unitary trans- 
form, e~ iKAt , gives an approximation to the time evolution on the 
lattice due to the kinetic energy operator that is accurate to ATth 
order in space. This leads to the digital quantum particle simulation 
algorithm: 

\xj,(t)) : = (e- iV *Fe- iTAt P*) ii \il/(p)). 



The QFT takes of order n 2 gates to calculate 55 and general algorithms 
implementing the diagonal T and V operators require ancillary 
qubits 56 , although it has been shown that the quadratic kinetic energy 
operator may be computed with n 2 two -qubit gates with no ancillary 
qubits 45 . 

The point of this paper is to describe a quantum particle simu- 
lation that can be used as a proof-of-principle demonstration. To do 
this, it must be implemented in current quantum computer archi- 
tectures where the number of qubits is limited and the number of 
gates that may be performed before decoherence destroys the com- 
putation is also limited. We consider the special case of square-well 
potentials. We show that a set of square-well potentials may be 
implemented with a sole, single-qubit operator and no ancillary 
qubits (see Methods). This virtually eliminates the calculation of 
the potential operator. Thus, simulations of important physical phe- 
nomena such as tunneling and the evolution of quasi-stable states 
may be performed with considerably fewer gates than simulations 
requiring arbitrary potentials. Additionally, we present a reduction 
of the number of gates required for the diagonal part of the kinetic 
operator in the case of 4 or fewer qubits. 

Results 

A Two Qubit Tunneling Simulation. Let us consider the smallest 
possible tunneling simulation. An n = 2 qubit simulation of N = 4 
lattice points may be performed with the circuit 



|Vo) ~ {#o} 



One Time Step of a Two-qubit Digital Quantum Single-particle 
Tunneling Simulation. 

In this circuit and below, Qy and Oy are controlled phase gates 
applied to qubits i and; and H f and Z f are a Hadamard gate and a Z- 
rotation on qubit i, respectively (see Methods). In total, each time 
step in this tunneling simulation requires 10 operations, 7 single 
qubit operations and 3 two-qubit operations. The circuit for a single 
time step is shown above. This circuit implements a double-well 
potential with the gate P 0 = exp( — iv<7 0 z At) acting on the lowest 
order qubit. 




In Fig. 1, we plot lattice occupation probabilities from two simula- 
tions with a double- well potential with At = 1/10: one a free-particle 
simulation with v = 0 and the other a tunneling simulation with v = 
10. This value of At traded off accuracy for gate number (see Methods 
section for a discussion of simulation errors), but the qualitative 
dynamics did not change even for more accurate (and costly) simu- 
lations. The initial state was |^/ m - ft -) = |01), corresponding to a par- 
ticle in one of the wells. The free-particle probability distribution 
spreads across all lattice points as it evolves, whereas the particle 
tunnels from the well at lattice point 1 (|01)) to the well at lattice 
point 3 (|11)) in the tunneling simulation. These results show that 
differences in the evolution of the probability distribution are evident 
within 4 time steps. Thus, such a proof-of-principle simulation may 
be implemented on a quantum computer with 4 X 10 = 40 gates (28 
single-qubit and 12 two-qubit). 

Multi-Qubit Quantum Tunneling Simulations. Larger simulations 
require more gates. For instance, a three-qubit simulation requires 6 
gates per QFT (3 singlequbit and 3 two-qubit), 6 gates for the 
diagonal kinetic energy operator (3 singlequbit and 3 two-qubit), 
and one single-qubit gate for the potential as shown in the circuit 
diagram below. 



I^o) 
l^i) 
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where Pi = exp ( — iva l z At) is shown, representatively, for a double- 
well potential, but other square-well potentials could be generated by 
acting on different qubits. The QFT and diagonal kinetic energy 
operators are 
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In Fig. 2, we show results from a three-qubit simulation with a 
double-well potential, where each well is resolved with two lattice 
points. The time step At = 1/5 and v = 5. Here, a tradeoff was also 
necessary to find a simulation that captured the oscillatory and tun- 
neling time scales and had few gates. The initial state was |i/^ m -#) = 
1 110). Because the initial state only occupies half of one well, oscil- 
latory dynamics are visible within the well. After a few time steps, the 
oscillatory state tunnels between wells. Oscillatory dynamics are 
evident within 4 time steps, but 5 time steps are required before 
the oscillatory state tunnels appreciably to the second well and 7 time 
steps are needed to see oscillation of the tunneled state. Thus, 
between 4 X 19 = 76 or 7X 19 = 133 gates would be required in 
order to see interesting tunneling effects in such a simulation. Four- 
qubit simulations can be envisaged using circuits based on similar 
methodology, although the number of gates would most likely be 
prohibitive for current quantum computer architectures. 

Discussion 

Only n = log 2 N qubits are required for an N lattice-point particle 
simulation, therefore this algorithm is efficient in the number of 
qubit resources required. In general, n 2 gates are needed for the 
kinetic operator (see Methods), but for fewer than n = 4 qubits, 



SCIENTIFIC REPORTS | 2 : 597 | DOI: 1 0.1 038/srep00597 



2 




Time 



10 V hi V l0 



Figure 2 | Particle probability distribution as a function of time over ten time steps in a three-qubit double-well simulation. The potential schematic is 
shown to the right, as in Fig. 1. Gray scales are as in Fig. 1. The double-well potential has two lattice points per well and, for the given initial state, an 
oscillation is induced in one of the wells then tunnels to the other well. 



further efficiency is possible using the basis approach that we out- 
lined. Furthermore, as noted above, our algorithm virtually elimi- 
nates the calculation of the potential operator. With very few qubits, 
interesting tunneling dynamics may be simulated with a gate count 
that is within reach of current quantum architectures. 

Much recent work has been done to understand the resources 
necessary to perform fully error- corrected quantum particle simula- 
tions using various error correction schemes 57 . We note that, as in the 
non-fault-tolerant case, the simulations presented here would also 
require the fewest resources in fully error corrected quantum archi- 
tectures, since the number of logical qubits and gates is as small as 
possible in the error-corrected case as well. 

Methods 

Square-well Potentials. A set of square-well potentials may be implemented with a 
sole single-qubit operator and no ancillary qubits. To see this, consider the single- 
qubit Z- rotation on the highest order qubit 

e -iVAt =e -iva" z - 1 At = e -ivo z At ( g )mi _ _ ? 

where v is a parameter, a superscript indicates the qubit to which the operator is 
applied and o z is the Pauli z- matrix a z =\ ® J .The operator e~ iVAt , when acting 



The kinetic operator is then K = FDF t , with 

D = exp ( - i( - 27i /4) 2 diag(0,4, 1 , 1 ) At) 



(1) 



on a lattice state, implements a square- well potential by rotating qubit states with |0) 
(|1> resp.) highest order qubit with positive (negative resp.) phase velocity v. The 
single-qubit operator acting on the next highest order qubit 



implements a double-square- well potential, and so on, with the last potential in this 
series implementing a Dirac comb-like potential. 

By simulating this class of square-well potentials we reduce the complexity of the 
potential calculation from a potentially large number of gates and ancillary qubits to 
just one single-qubit operation. In contrast, the phase kickback algorithm used in 41 
requires a total of 2n qubits (n for the quantum state and n ancilla qubits) as well as n 2 
extra gates for a QFT on n ancilla qubits. 

Although the main point of this paper is the simplification of the simulation by 
reducing the potential computation to one gate, further gate reductions can be made 
to the kinetic operator when considering few qubit systems. This can be done by 
forming a diagonal Hamiltonian using a basis of diagonal operators. This method is 
only efficient for n < 4 and scales as 2" — 1. Therefore, for n > 4, ancillary qubits or 
the Benenti-Strini (BS) 45 method (that scales as n 2 ) should be used. However, for small 
systems of qubits, forming a basis of diagonal operators takes fewer gates. Note that 
for n = 2, 3 and 4, we would need 4, 9 and 16 (resp.) operators for BS, but at most 3, 7 
and 15 (resp.) are required using a basis along the diagonal. Furthermore, fewer multi- 
qubit operations are necessary with this method, although a three-qubit operator may 
be necessary for three-qubit simulations and three- and four-qubit operators for four- 
qubit simulations. Here, we use the method to construct the diagonal operator D for 
the kinetic energy operator, however, arbitrary potentials could also be constructed 
this way. Because our goal is to find interesting simulations with few gates, we do not 
pursue more complex potentials here. 

In our two-qubit simulations (see Results), we use a double well potential, 
exp( — iVAt) = exp ( — iva° z At) . The QFT may be computed with the operators 55 

F^ =HiQ 0 iHq, 

where is a Hadamard operator on the z'th qubit and the controlled-phase gate Q 0 i 
= diag(l, 1, 1, co) with co = exp(27n'/4). Note that F^ results in a bitswapped Fourier 
transform (in our notation, F is the inverse Fourier transform matrix). 



(Note that this operator is also bit- swapped and we have taken m = 1/2). 

The diagonal operation D may be achieved up to an overall phase with the operator 



D = OoiZiZo 



where the single-qubit operators 



and O 01 is a controlled-phase operator on qubits 0 and 1, 



iyc 2 diag(l,l,l,-l) 01 Af 



where y = ( — 2ti/4) 2 / y/4. The coefficients (c 0 = — 1, q = —4 and c 2 = 4) in the unitary 
operators Z 0 , Z x and O 01 (resp.) were obtained by noting that the vectors (1, 1,— 1,— 1), 
(1,-1,1,-1) and (1, 1, 1, — 1) that form their diagonal elements are a basis for zero mean 
vectors (i.e. neglecting the constant phase proportional to (1, 1, 1, 1)) for U 4 . 

To calculate the QFT in the three-qubit case, the unitary operators are the single-qubit 
Hadamard gates and the two-qubit controlled-phase operators Qy = Qy(co), where co = 
exp(27n'/8) (see circuit diagram in Results). To calculate the quadratic diagonal operator 
in the three-qubit case (see Results), the unitary operators are the single-qubit operators 

Zo = exp(-iyc 0 (j 0 z At) 

Zi = exp ( — iyci o\At) 

Z 2 =exp( — iyc 2 CT 2 z At) 
and the two-qubit controlled-phase operators 

001 = exp ( - iyc 3 diag( 1 , 1, 1 , - 1) 01 At) 

0 02 = exp ( - iy c 4 diag( 1 , 1 , 1 , - 1 ) ^ At) 
0 12 =exp(-zyc 5 diag(l,l,l,-l) 12 Af), 

where 

y=-(27i/8)Yv / 8 

and 

Co =-1.42 
Cl = -5.66 
c 2 = - 22.63 

c 3 = 22.63 

c 4 = 11.31 
c 5 = -5.66. 

Note that, in principle, a seventh three-qubit operator would also be necessary propor- 
tional to the vector (1, 1, 1, 1, 1, 1, 1,-1) along the diagonal, but for the (bit- swapped) 
diagonal of the 8 lattice point simulation, (0, 16, 4, 4, 1, 9, 9, 1), its coefficient is identically 
zero. This circuit requires a total of 19 gates (10 single-qubit and 9 two-qubit). 

Once we have computed the unitary operators for one time step of the kinetic, 
U k (At), and potential, U v (At), evolution, it is straightforward to calculate the exact 
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Table 1 | RMS Error per time step for a range of values of At for the 
2 and 3 qubit simulations described in Results 

At 0.2 0.1 0.05 0.025 0.0125 



Error 2 qubit 0.40 0.16 0.06 0.03 0.01 
Error 3 qubit 0.08 0.02 0.007 0.003 0.002 



Hamiltonian evolution: J7 exact (Af) = expflog([7*(Af))+log(l7 v (Af))). We then 
determine RMS errors relative to the exact evolution for the two and three-qubit 
simulations. These are given in Table 1. The simulations that we present in Figs. 1 and 
2 have a total, integrated error of 64% and 20%, respectively. For both simulations, we 
checked that state evolution was qualitatively similar to the exact simulations for the 
chosen values of At. Clearly, for high-precision simulations, smaller errors would be 
desirable. To achieve this, one would use smaller values for At at the expense of an 
increase in gate count. 
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